	SUBROUTINE YUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA)
	REAL(8),DIMENSION(L,M,N)::U,V,DIV,OMEGA,OMEGA1
	REAL(8),DIMENSION(L,M,N)::OMT	 !OMT为由热力学法算出的顶层OMEGA值
	REAL(8),DIMENSION(N)::P
	REAL(8)::F0,DL,DM
	OMT=0   !这里取OMT=0,如有资料可用热力学法计算
	OMEGA1=0
	CALL DIVER(U,V,P,DL,DM,F0,L,M,N,DIV)
	DO I=1,L
	  DO J=1,M
	    DO K=2,N
	OMEGA1(I,J,K)=OMEGA1(I,J,K-1)-(DIV(I,J,K-1)+DIV(I,J,K))/2*(P(K)-P(K-1))
		END DO
	    DO K=1,N
	OMEGA(I,J,K)=OMEGA1(I,J,K)-K*(K-1)/(N*(N-1))*(OMEGA1(I,J,N)-OMT(I,J,N))
		END DO
	  END DO
	END DO
	END

	SUBROUTINE OME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA)  !求解准地转OMEGA方程
	REAL(8),DIMENSION(L,M,N)::U,V,FAI,KAP,SIGMA
	REAL(8),DIMENSION(L,M,5)::T,TP
	REAL(8),DIMENSION(L,M,N)::OMEGA1,OMEGA2,OMEGA,OMEGA0
	REAL(8),DIMENSION(L,M,N)::QA1,QA2,QA3,QB1,QB2,QB3	 !中间变量
	REAL(8),DIMENSION(N)::P,SIGMAV
	REAL(8),DIMENSION(M)::FCO
	REAL(8)::EPS,CP,RD,F0,DL,DM,ALF
	EPS=1.E-5   !迭代要求的精度
	ALF=1.44
	RD=287
	CP=1005
	OMEGA1=0
	OMEGA2=0
	CALL LHS(U,V,P,L,M,N,F0,DL,DM,FAI,KAP)
	DO J=1,M		!科氏参数
	FCO(J)=2*(7.29E-5)*SIN(F0+(J-1)*DM)
	END DO
!   计算静力稳定度
	CALL AP1(T,L,M,N,P,TP)
	DO K=1,N
	  DO J=1,M
	    DO I=1,L
	      SIGMA(I,J,K)=-(RD/P(K)*(TP(I,J,K)-RD/CP*T(I,J,K)/P(K)))
	    END DO
	  END DO
	END DO
	DO K=1,N
      SIGMAV(K)=0
	  DO I=1,L
	   DO J=1,M
		 SIGMAV(K)=SIGMAV(K)+SIGMA(I,J,K)
	   END DO
     	END DO
   	END DO
	DO K=1,N
	  DO J=1,M
	    DO I=1,L
 	      IF(SIGMA(I,J,K)<=0)THEN
			SIGMA(I,J,K)=SIGMAV(K)
	      END IF
	    END DO
      END DO
	END DO
!	计算右边第一项	
	CALL LAPLACE(FAI,L,M,N,F0,DL,DM,QA1)
	DO K=1,N
	  DO J=1,M
	    DO I=1,L
	      QA1(I,J,K)=QA1(I,J,K)+FCO(J)
	    END DO
	  END DO
	END DO
	CALL JACOB(FAI,QA1,DL,DM,F0,L,M,N,QA2)
	CALL AP1(QA2,L,M,N,P,QA3)
	QA3=QA3*FCO(M/2)
	N0=0
	DO 
	  N0=N0+1
	  OMEGA0=OMEGA1
	  CALL ERROR(OMEGA1,QA3,P,SIGMA,L,M,N,F0,FCO,DL,DM,ALF)
	  IF(MAXVAL(ABS(OMEGA1-OMEGA0))<EPS)EXIT
	END DO
	open(2,file='omega1.dat')
	WRITE(2,'(<L>E11.4)')(((OMEGA1(I,J,K),I=1,L),J=1,M),K=1,N)
	close(2)
!	计算右边第二项	
	CALL AP1(FAI,L,M,N,P,QB1)
	CALL JACOB(FAI,QB1,DL,DM,F0,L,M,N,QB2)
	CALL LAPLACE(QB2,L,M,N,F0,DL,DM,QB3)
	QB3=-FCO(M/2)*QB3
	N1=0
	DO 
	  N1=N1+1
	  OMEGA0=OMEGA2
	  CALL ERROR(OMEGA2,QB3,P,SIGMA,L,M,N,F0,FCO,DL,DM,ALF)
	  IF(MAXVAL(ABS(OMEGA2-OMEGA0))<EPS)EXIT
	END DO
	open(3,file='omega2.dat')
	WRITE(3,'(<L>E11.4)')(((OMEGA2(I,J,K),I=1,L),J=1,M),K=1,N)
	close(3)
	OMEGA=OMEGA1+OMEGA2
	END 

	SUBROUTINE AP1(A,L,M,N,P,FF)		 !计算垂直方向一次不等距差分
	REAL(8),DIMENSION(L,M,N)::A
	REAL(8),DIMENSION(L,M,N)::FF
	REAL(8),DIMENSION(N)::P
	DO J=1,M
	  DO I=1,L
	    DO K=2,N-1
	FF(I,J,K)=((P(K)-P(K-1))*A(I,J,K+1)+(P(K+1)+P(K-1)-2*P(K))*A(I,J,K)
     &-(P(K+1)-P(K))*A(I,J,K-1))/(2*(P(K)-P(K-1))*(P(K+1)-P(K)))
        END DO	
	  END DO	
    	END DO	
	END

	SUBROUTINE LAPLACE(A,L,M,N,F0,DL,DM,QA)
	REAL(8),DIMENSION(L,M,N)::A,QA
	REAL(8)::F0,DL,DM,AA
	AA=6.371E6
	L1=L-1
	M1=M-1
	DO K=1,N
	  DO I=2,L1
	    DO J=2,M1
		  QA(I,J,K)=(A(I+1,J,K)+A(I-1,J,K)-2*A(I,J,K))/(DL*COS(F0+(J-1)
     &*DM))**2+(A(I,J+1,K)+A(I,J-1,K)-2*A(I,J,K))/(DM*DM)
     &-TAN(F0+(J-1)*DM)*(A(I,J+1,K)-A(I,J-1,K))/(2*DM)
		END DO
	  END DO
	END DO
	QA=QA/(AA*AA)
	CALL BOUND(QA,L,M,N)
	END

	SUBROUTINE JACOB(A,B,DL,DM,F0,L,M,N,JC)  !采用两种JACOB差分方法的平均
	REAL(8),DIMENSION(L,M,N)::A,B,JC
	REAL(8)::F0,DL,DM,J1,J2     
	REAL(8)::AA=6.371E6
	L1=L-1
	M1=M-1
	DO K=1,N
	  DO J=2,M1
	    DO I=2,L1
			J1=-((B(I+1,J,K)-B(I-1,J,K))*(A(I,J+1,K)-A(I,J-1,K))
     &-(B(I,J+1,K)-B(I,J-1,K))*(A(I+1,J,K)-A(I-1,J,K)))
     &/(AA*AA*COS(F0+(J-1)*DM)*4*DL*DM)
			J2=-(B(I+1,J+1,K)*(A(I,J+1,K)-A(I+1,J,K))
     &-B(I-1,J-1,K)*(A(I-1,J,K)-A(I,J-1,K))
     &-B(I-1,J,K)*(A(I,J+1,K)-A(I-1,J,K))
     &+B(I+1,J-1,K)*(A(I+1,J,K)-A(I,J-1,K)))
     &/(AA*AA*COS(F0+(J-1)*DM)*4*DL*DM)
		  JC(I,J,K)=0.6*J1+0.4*J2
		END DO
	  END DO
	END DO
	CALL BOUND(JC,L,M,N)
	END

	SUBROUTINE ERROR(A,Q,P,SI,L,M,N,F0,FCO,DL,DM,ALF)
	REAL(8),DIMENSION(L,M,N)::A,Q,SI
	REAL(8),DIMENSION(N)::P
	REAL(8),DIMENSION(M)::FCO
	REAL(8)::F0,DL,DM,AA,EPS,ALF,R
	EPS=1.E3
	AA=6.371E6
	L1=L-1
	M1=M-1
	N1=N-1
	DO K=2,N1
	 DO I=2,L1
	   DO J=2,M1
		   R=((A(I+1,J,K)+A(I-1,J,K))/(DL*COS(F0+(J-1)*DM))**2
     &	   +(A(I,J+1,K)+A(I,J-1,K))/DM**2-TAN(F0+(J-1)*DM)*
     &		  (A(I,J+1,K)-A(I,J-1,K))/(2*DM)+(AA*FCO(M/2))**2*
     &		  ((P(K)-P(K-1))*A(I,J,K+1)+(P(K+1)-P(K))*A(I,J,K-1))
     &		  /(SI(I,J,K)*(P(K)-P(K-1))**2*(P(K+1)-P(K)))
     &		  -AA*AA*Q(I,J,K)/SI(I,J,K))/(2/(DL*COS(F0+(J-1)*DM))**2
     &		  +2/DM**2+(AA*FCO(M/2))**2*(P(K+1)-P(K-1))
     &		  /(SI(I,J,K)*(P(K)-P(K-1))**2*(P(K+1)-P(K))))
	   A(I,J,K)=ALF*R+(1-ALF)*A(I,J,K)
	  END DO
	 END DO
	END DO
	END
	
	PROGRAM MAIN
	INTEGER,PARAMETER::L=57,M=29,N0=7,N=10
	REAL(8),DIMENSION(L,M,N0)::U0,V0,T0
	REAL(8),DIMENSION(N0)::P0
	REAL(8),DIMENSION(L,M,N)::U,V,T
	REAL(8),DIMENSION(L,M,N)::OMEGA
	REAL(8),DIMENSION(N)::P
	REAL(8)::F0,DL,DM,PI
	DATA P0/100,200,300,500,700,850,1000/
	DATA P/1000,900,800,700,600,500,400,300,200,100/
	L1=L
	PI=3.1415926
	F0=0*PI/180.
	DL=2.5*PI/180
	DM=2.5*PI/180
	WRITE(*,'(" INPUT IO(IO=1,用运动学方法；IO=2,用OMEGA方程)")')
	READ(*,'(I1)')IO
	OPEN(12,FILE='U880501.DAT')
	READ(12,'(<L1>F7.1)')(((U0(I,J,K),I=1,L),J=1,M),K=N0,1,-1)
	CLOSE(12)
	OPEN(13,FILE='V880501.DAT')
	READ(13,'(<L1>F7.1)')(((V0(I,J,K),I=1,L),J=1,M),K=N0,1,-1)
	CLOSE(13)
	OPEN(13,FILE='T880501.DAT')
	READ(13,'(<L1>F7.1)')(((T0(I,J,K),I=1,L),J=1,M),K=N0,1,-1)
	CLOSE(13)
	T0=T0+273.15
!	垂直插值
	DO I=1,L
	  DO J=1,M
	    CALL SPLINECALCULATE(P0,U0(I,J,1:N0),N0)
	    DO K=1,N
		  CALL SPLINEEVALUATE(P(K),U(I,J,K))
	    END DO	
	    CALL SPLINECALCULATE(P0,V0(I,J,1:N0),N0)
	    DO K=1,N
		  CALL SPLINEEVALUATE(P(K),V(I,J,K))
	    END DO	
	    CALL SPLINECALCULATE(P0,T0(I,J,1:N0),N0)
	    DO K=1,N
		  CALL SPLINEEVALUATE(P(K),T(I,J,K))
	    END DO	
	  END DO
	END DO
	IF(IO==1)THEN
	CALL YUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA)
	OPEN(20,FILE='OMEGAY.DAT')
	WRITE(20,'(<L1>E11.4)')(((OMEGA(I,J,K),I=1,L),J=1,M),K=1,N)
	CLOSE(20)
	ELSE
	CALL OME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA)
	OPEN(20,FILE='OMEGA.DAT')
	WRITE(20,'(<L1>E11.4)')(((OMEGA(I,J,K),I=1,L),J=1,M),K=1,N)
	CLOSE(20)
	ENDIF
	END

! SPLINE.FOR -- NATURAL CUBIC SPLINE
	SUBROUTINE SPLINECALCULATE(INX,INA,NUMX )
	IMPLICIT NONE
	REAL(8)::INA(*),INX(*)
	INTEGER::NUMX
	INTEGER::MAXPOINTS
	PARAMETER(MAXPOINTS=1000 )
	INTEGER::NXS
	REAL(8)::X(0:MAXPOINTS)
	REAL(8)::A(0:MAXPOINTS),B(0:MAXPOINTS)
	REAL(8)::C(0:MAXPOINTS),D(0:MAXPOINTS)
	COMMON /SPLINEDATA/ A, B, C, D, NXS, X
	INTEGER::I
	REAL(8)::ALPHA(0:MAXPOINTS), L(0:MAXPOINTS)
	REAL(8)::MU(0:MAXPOINTS), Z(0:MAXPOINTS)
	REAL(8)::H(0:MAXPOINTS)
	NXS=NUMX-1	
	DO I=0, NXS
		X(I)=INX(I+1)
		A(I)=INA(I+1)
	END DO
	DO I=0, NXS-1
		H(I)=X(I+1)-X(I)
	END DO
	X(NXS+1)=40000
	DO I=1, NXS-1
		ALPHA(I)=3.0*(A(I+1)*H(I-1)-A(I)*
     &		     (X(I+1)-X(I-1))+A(I-1)*H(I))/(H(I-1)*H(I))
	END DO
	L(0)=1.0
	MU(0)=0.0
	Z(0)=0.0
	DO I=1, NXS-1
		L(I)=2.0*(X(I+1)-X(I-1))-H(I-1)*MU(I-1)
		MU(I)=H(I)/L(I)
		Z(I)=(ALPHA(I)-H(I-1)*Z(I-1))/L(I)	
	END DO
	L(NXS)=1.0
	Z(NXS)=0.0
	C(NXS)=0.0
	DO I=NXS-1, 0, -1
		C(I)=Z(I)-MU(I)*C(I+1)
		B(I)=(A(I+1)-A(I))/H(I)-H(I)*(C(I+1)+2.0*C(I))/3.0
		D(I)=(C(I+1)-C(I))/(3.0*H(I))
	END DO
	END

	SUBROUTINE SPLINEEVALUATE(EVALX,EVALY )	
	REAL(8)::EVALX, EVALY
	INTEGER:: MAXPOINTS
	PARAMETER( MAXPOINTS=1000 )
	INTEGER::NXS
	REAL(8)::X(0:MAXPOINTS)
	REAL(8)::A(0:MAXPOINTS), B(0:MAXPOINTS)
	REAL(8)::C(0:MAXPOINTS), D(0:MAXPOINTS)
	COMMON /SPLINEDATA/ A, B, C, D, NXS, X
	REAL(8)::TERM
	INTEGER::I
	I=0	
	DO WHILE(.NOT.((X(I).LE.EVALX).AND.(X(I+1).GT.EVALX)))
		I=I+1
	END DO	
	TERM=EVALX-X(I)
        EVALY=A(I)+B(I)*TERM+C(I)*TERM*TERM+D(I)*TERM*TERM*TERM
	END
